Multicontrast MRI reconstruction with structure-guided total variation

Magnetic resonance imaging (MRI) is a versatile imaging technique that allows different contrasts depending on the acquisition parameters. Many clinical imaging studies acquire MRI data for more than one of these contrasts---such as for instance T1 and T2 weighted images---which makes the overall scanning procedure very time consuming. As all of these images show the same underlying anatomy one can try to omit unnecessary measurements by taking the similarity into account during reconstruction. We will discuss two modifications of total variation---based on i) location and ii) direction---that take structural a priori knowledge into account and reduce to total variation in the degenerate case when no structural knowledge is available. We solve the resulting convex minimization problem with the alternating direction method of multipliers that separates the forward operator from the prior. For both priors the corresponding proximal operator can be implemented as an extension of the fast gradient projection method on the dual problem for total variation. We tested the priors on six data sets that are based on phantoms and real MRI images. In all test cases exploiting the structural information from the other contrast yields better results than separate reconstruction with total variation in terms of standard metrics like peak signal-to-noise ratio and structural similarity index. Furthermore, we found that exploiting the two dimensional directional information results in images with well defined edges, superior to those reconstructed solely using a priori information about the edge location.Engineering and Physical Sciences Research Council (Grant ID: EP/H046410/1)This is the final version of the article. It first appeared from Society for Industrial and Applied Mathematics via http://dx.doi.org/10.1137/15M1047325

Mathematics of biomedical imaging today—a perspective

Biomedical imaging is a fascinating, rich and dynamic research area, which has huge importance in biomedical research and clinical practice alike. The key technology behind the processing, and automated analysis and quantification of imaging data is mathematics. Starting with the optimisation of the image acquisition and the reconstruction of an image from indirect tomographic measurement data, all the way to the automated segmentation of tumours in medical images and the design of optimal treatment plans based on image biomarkers, mathematics appears in all of these in different flavours. Non-smooth optimisation in the context of sparsity-promoting image priors, partial differential equations for image registration and motion estimation, and deep neural networks for image segmentation, to name just a few. In this article, we present and review mathematical topics that arise within the whole biomedical imaging pipeline, from tomographic measurements to clinical support tools, and highlight some modern topics and open problems. The article is addressed to both biomedical researchers who want to get a taste of where mathematics arises in biomedical imaging as well as mathematicians who are interested in what mathematical challenges biomedical imaging research entails

Hamilton-Green solver for the forward and adjoint problems in photoacoustic tomography

The majority of the solvers for the acoustic problem in Photoacoustic Tomography (PAT) rely on full solution of the wave equation, which makes them less suitable for real-time and dynamic applications where only partial data is available. This is in contrast to other tomographic modalities, e.g. X-ray tomography, where partial data implies partial cost for the application of the forward and adjoint operators. In this work we present a novel solver for the forward and adjoint wave equations for the acoustic problem in PAT. We term the proposed solver Hamilton-Green as it approximates the fundamental solution to the respective wave equation along the trajectories of the Hamiltonian system resulting from the high frequency asymptotic approximate solution for the wave equation. This approach is flexible and scalable in the sense that it allows computing the solution for each sensor independently at a fraction of the cost of the full wave solution. The theoretical foundations of our approach are rooted in results available in seismics and ocean acoustics. To demonstrate the feasibility of our approach we present results for 2D domains with homogeneous and heterogeneous sound speeds and evaluate them against a full wave solution obtained with a pseudospectral finite difference method implemented in the k-Wave toolbox [1]

Restarting iterative projection methods for Hermitian nonlinear eigenvalue problems with minmax property

In this work we present a new restart technique for iterative projection methods for nonlinear eigenvalue problems admitting minmax characterization of their eigenvalues. Our technique makes use of the minmax induced local enumeration of the eigenvalues in the inner iteration. In contrast to global numbering which requires including all the previously computed eigenvectors in the search subspace, the proposed local numbering only requires a presence of one eigenvector in the search subspace. This effectively eliminates the search subspace growth and therewith the super-linear increase of the computational costs if a large number of eigenvalues or eigenvalues in the interior of the spectrum are to be computed. The new restart technique is integrated into nonlinear iterative projection methods like the Nonlinear Arnoldi and Jacobi-Davidson methods. The efficiency of our new restart framework is demonstrated on a range of nonlinear eigenvalue problems: quadratic, rational and exponential including an industrial real-life conservative gyroscopic eigenvalue problem modeling free vibrations of a rolling tire. We also present an extension of the method to problems without minmax property but with eigenvalues which have a dominant either real or imaginary part and test it on two quadratic eigenvalue problems

Multi-sheet surface rebinning methods for reconstruction from asymmetrically truncated cone beam projections: I. Approximation and optimality

The mechanical motion of the gantry in conventional cone beam CT scanners restricts the speed of data acquisition in applications with near real time requirements. A possible resolution of this problem is to replace the moving source detector assembly with static parts that are electronically activated. An example of such a system is the Rapiscan Systems RTT80 real time tomography scanner, with a static ring of sources and axially offset static cylinder of detectors. A consequence of such a design is asymmetrical axial truncation of the cone beam projections resulting, in the sense of integral geometry, in severely incomplete data. In particular we collect data only in a fraction of the Tam-Danielsson window, hence the standard cone beam reconstruction techniques do not apply. In this work we propose a family of multi-sheet surface rebinning methods for reconstruction from such truncated projections. The proposed methods combine analytical and numerical ideas utilizing linearity of the ray transform to reconstruct data on multi-sheet surfaces, from which the volumetric image is obtained through deconvolution. In this first paper in the series, we discuss the rebinning to multi-sheet surfaces. In particular we concentrate on the underlying transforms on multi-sheet surfaces and their approximation with data collected by offset multi-source scanning geometries like the RTT. The optimal multi-sheet surface and the corresponding rebinning function are found as a solution of a variational problem. In the case of the quadratic objective, the variational problem for the optimal rebinning pair can be solved by a globally convergent iteration. Examples of optimal rebinning pairs are computed for different trajectories. We formulate the axial deconvolution problem for the recovery of the volumetric image from the reconstructions on multi-sheet surfaces. Efficient and stable solution of the deconvolution problem is the subject of the second paper in this series (Betcke and Lionheart 2013 Inverse Problems 29 115004). © 2013 IOP Publishing Ltd

Iterated preconditioned LSQR method for inverse problems on unstructured grids

This article presents a method for solving large-scale linear inverse imaging problems regularized with a nonlinear, edge-preserving penalty term such as total variation or the Perona–Malik technique. Our method is aimed at problems defined on unstructured meshes, where such regularizers naturally arise in unfactorized form as a stiffness matrix of an anisotropic diffusion operator and factorization is prohibitively expensive. In the proposed scheme, the nonlinearity is handled with lagged diffusivity fixed point iteration, which involves solving a large-scale linear least squares problem in each iteration. Because the convergence of Krylov methods for problems with discontinuities is notoriously slow, we propose to accelerate it by means of priorconditioning (Bayesian preconditioning). priorconditioning is a technique that, through transformation to the standard form, embeds the information contained in the prior (Bayesian interpretation of a regularizer) directly into the forward operator and thence into the solution space. We derive a factorization-free preconditioned LSQR algorithm (MLSQR), allowing implicit application of the preconditioner through efficient schemes such as multigrid. The resulting method is also matrix-free i.e. the forward map can be defined through its action on a vector. We illustrate the performance of the method on two numerical examples. Simple 1D-deblurring problem serves to visualize the discussion throughout the paper. The effectiveness of the proposed numerical scheme is demonstrated on a three-dimensional problem in fluorescence diffuse optical tomography with total variation regularization derived algebraic multigrid preconditioner, which is the type of large scale, unstructured mesh problem, requiring matrix-free and factorization-free approaches that motivated the work here

Photoacoustic Reconstruction Using Sparsity in Curvelet Frame: Image Versus Data Domain

Curvelet frame is of special significance for photoacoustic tomography (PAT) due to its sparsifying and microlocalisation properties. We derive a one-to-one map between wavefront directions in image and data spaces in PAT which suggests near equivalence between the recovery of the initial pressure and PAT data from compressed/subsampled measurements when assuming sparsity in Curvelet frame. As the latter is computationally more tractable, investigation to which extent this equivalence holds conducted in this paper is of immediate practical significance. To this end we formulate and compare DR , a two step approach based on the recovery of the complete volume of the photoacoustic data from the subsampled data followed by the acoustic inversion, and p0R , a one step approach where the photoacoustic image (the initial pressure, p0 ) is directly recovered from the subsampled data. Effective representation of the photoacoustic data requires basis defined on the range of the photoacoustic forward operator. To this end we propose a novel wedge-restriction of Curvelet transform which enables us to construct such basis. Both recovery problems are formulated in a variational framework. As the Curvelet frame is heavily overdetermined, we use reweighted ℓ1 norm penalties to enhance the sparsity of the solution. The data reconstruction problem DR is a standard compressed sensing recovery problem, which we solve using an ADMM-type algorithm, SALSA. Subsequently, the initial pressure is recovered using time reversal as implemented in the k-Wave Toolbox. The p0 reconstruction problem, p0R , aims to recover the photoacoustic image directly via FISTA, or ADMM when in addition including a non-negativity constraint. We compare and discuss the relative merits of the two approaches and illustrate them on 2D simulated and 3D real data in a fair and rigorous manner

Acoustic Wave Field Reconstruction From Compressed Measurements With Application in Photoacoustic Tomography

We present a method for the recovery of compressively sensed acoustic fields using patterned, instead of point-by-point, detection. From a limited number of such compressed measurements, we propose to reconstruct the field on the sensor plane in each time step independently assuming its sparsity in a Curvelet frame. A modification of the Curvelet frame is proposed to account for the smoothing effects of data acquisition and motivated by a frequency domain model for photoacoustic tomography. An ADMM type algorithm, split augmented Lagrangian shrinkage algorithm, is used to recover the pointwise data in each individual time step from the patterned measurements. For photoacoustic applications, the photoacoustic image of the initial pressure is reconstructed using time reversal in k-Wave Toolbox

Feasibility Study of Time-of-Flight Compton Scatter Imaging Using Picosecond Length X-Ray Pulses

By measuring the time of flight of scattered X-ray photons, the point of interaction, assuming a single scatter, can be determined, providing 3-D information about an object under inspection. This paper describes experimental ToF Compton scatter measurements conducted at the versatile electron linear accelerator (VELA), a picosecond pulsewidth electron source situated in Daresbury, U.K. The ToF of scattered X-ray photons was measured using a CeBr3 detector, and a full width at half maximum (FWHM) of between 29 and 36 cm was achieved with a 5-cm-thick plastic test object. By implementing a low-energy cutoff, the FWHM was reduced to between 12 and 26 cm. Two test objects placed in series with a 50-cm space between were separable in the data after applying the low energy cutoff

Mathematics of biomedical imaging today—a perspective

Abstract Biomedical imaging is a fascinating, rich and dynamic research area, which has huge importance in biomedical research and clinical practice alike. The key technology behind the processing, and automated analysis and quantification of imaging data is mathematics. Starting with the optimisation of the image acquisition and the reconstruction of an image from indirect tomographic measurement data, all the way to the automated segmentation of tumours in medical images and the design of optimal treatment plans based on image biomarkers, mathematics appears in all of these in different flavours. Non-smooth optimisation in the context of sparsity-promoting image priors, partial differential equations for image registration and motion estimation, and deep neural networks for image segmentation, to name just a few. In this article, we present and review mathematical topics that arise within the whole biomedical imaging pipeline, from tomographic measurements to clinical support tools, and highlight some modern topics and open problems. The article is addressed to both biomedical researchers who want to get a taste of where mathematics arises in biomedical imaging as well as mathematicians who are interested in what mathematical challenges biomedical imaging research entails.</jats:p